!===========================================================
! Test code to investigate how to accelerate the convergence
! of iterative methods for the solution of the Poisson 
! equation. Two-dimensional case
! The geometry is a Shelf-Break type, with a shallow and 
! a deep end.
! Uses routines from ELLIKIT. 
! Copyright 2008 A. Scotti
! Distributed under the GPL.
!=========================================================== 

USE ElliDef ! defines the data-type associated with the solvers in ELLIKIT
USE NodeInfoDef ! defines the data-type NodeInfo
USE metric ! contains the elements of the metric coefficient
USE Elliptic_Solvers ! Contains the main routines for generation and 

                     ! solution of elliptic problems
IMPLICIT NONE ! To avoid Snafus
TYPE (NodeInfo),DIMENSION(2) :: Info ! Collects all relevant infos.
                                     ! strictly speaking not necessary
                                     ! ELLIKIT only sees Info%ElliInfo


INTEGER error,N ! used to collect error status and size of problem
REAL (KIND=qPrec) :: xUpper(1:3),dx,dy,dz ! define computational space
REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: p,derp,u,v,w,div,pl ! sundry local variables
                                                                     
REAL (KIND=qPrec), ALLOCATABLE, DIMENSION (:,:,:) :: ut,vt,wt ! ditto


INTEGER i,j,k,nx,ny,nz,l ! local integers for indexes and size of problem
INTEGER, PARAMETER ::  how_many_steps=500 ! frequency of acceleration
! solver is used to set the particular solver for the iterative part. Only uncomment one at a time
!INTEGER, PARAMETER :: solver=PARDISO_SOLVER ! this is the direct solver. Not needed here
!INTEGER, PARAMETER :: solver=IT_SD ! conjugate gradient
!INTEGER, PARAMETER :: solver=IT_CG ! conjugate gradient
!INTEGER, PARAMETER :: solver=IT_BCG ! do not use, since the problem 
                                     ! is symmetric
!INTEGER, PARAMETER :: solver=IT_DBCG ! same as bcg, 
!INTEGER, PARAMETER :: solver=IT_BCGSTAB !BiConjugate Gradient stabilized
!INTEGER, PARAMETER :: solver=IT_TFQMR !TransposeFree QuasiMinimim Residual
INTEGER, PARAMETER :: solver=IT_FOM ! Fully Orthogonal
!INTEGER, PARAMETER :: solver=IT_GMRES ! Generalized Minimum Residual method
!INTEGER, PARAMETER :: solver=IT_FGMRES ! Flexible version of GMRES
!INTEGER, PARAMETER :: solver=IT_DQGMRES ! Direct version of Quasi GMRES

!INTEGER, PARAMETER :: solver=DCG_SOLVER ! same as CG but uses the Intel MKL version. Preconditioned version sucks
LOGICAL iterate ! used to decide when to stop.


! 0. Set up the problem
! Notes: 
! 1) we use two info structures. One to go with the iterative method, the other 
! to be used for the expansion method. The two are linked (see below) so that they share
! the same vector field. In real life, one structure should be enough
! 2) What really needs to be initialized here is Info%ElliInfo%... 
! The other parameters are really never used...

! populate Info structures
DO i=1,2 ! two Infos
   ALLOCATE(Info(i)%ElliInfo) ! ElliInfo is defined as a pointer, space needs to be allocated

   Info(i)%mX=(/50,140,1/) ! size of the numerical grid (nx,ny,nz). nz=1 since we are dealing with a two-D problem
   nx=Info(i)%mX(1);ny=Info(i)%mX(2);nz=Info(i)%mX(3) 
   Info(i)%nDim=2 ! two-D problem
   Info(i)%ElliInfo%ndim=Info(i)%ndim ! transfer the info to ElliInfo
   Info(i)%ElliInfo%mX=Info(i)%mX ! last dim has to be one
   
   Info(i)%xLower=(/ 1.d4/2.2d3,0.d0,0.d0 /) ! lower boundaries of logical domain
   
   Info(i)%xUpper=(/15000./2200,1.,1./) ! upper boundaries of logical domain
   
   Info(i)%dX(1:2)=(Info(i)%xUpper(1:2)-Info(i)%xLower(1:2))/Info(i)%mX(1:2);Info(i)%dX(3)=1 ! calculates dx,dy,dz
   Info(i)%ElliInfo%dX(1:2)=Info(i)%dX(1:2) ! transfer to ElliInfo
   Info(i)%ElliInfo%dX(3)=Info(i)%dX(3) ! just to have it set it. 
   Info(i)%dmax=0.
   
   ! assume all Neumann conditions
   Info(i)%ElliInfo%mthbc(1:6)=0 ! it means that the projection is done assigning normal velocities at boundaries.

   
   
   
   
   CALL setprob(Info(i)) ! Initializes the metric 
   PRINT*,'Dimension of problem',Info(i)%mX
   
   
   CALL DEFINE_TOPO(Info(i)) ! Calculates the metric coefficients (see metric.f90)
   PRINT*,'Max Depth =',MAXVAL(1-Info(i)%topo%h)*Info%dmax
   IF(i==1) THEN ! first Info is used to solve iterative problem
      error=CREATE_OPERATORS(Info(i)%ElliInfo)
   ELSE
      error=CREATE_OPERATORS_THIN(Info(i)%ElliInfo) ! second to implement the accelerator
   END IF
   ! zero BC here to be safe 
   DO j=1,2*Info(i)%nDim
      Info(i)%ElliInfo%bc(j)%p=zero ! stores bc
   END DO
   
   IF(.NOT.check(Info(i)%ElliInfo,CHECK_BC)) THEN ! this is to make sure that everything is ok. Redundant
      PRINT*,'BC are not consistent'
      STOP
   END IF
   
END DO ! end of initialization

! 1. Now we fill up the velocity fields with random numbers



ALLOCATE(ut(1:ny,1:nz,2),vt(1:nx,1:nz,2)) ! these are the values of the velocity at boundary
CALL RANDOM_NUMBER(ut) ! fill them with random numbers
CALL RANDOM_NUMBER(vt)
ut(:,1,1)=ut(:,1,1)-SUM(ut(:,1,1))/ny ! now make sure that flux across boundary is conserved 
ut(:,1,2)=ut(:,1,2)-SUM(ut(:,1,2))/ny ! or else all hell breaks loose.
vt(:,1,1)=vt(:,1,1)-SUM(vt(:,1,1))/nx
vt(:,1,2)=vt(:,1,2)-SUM(vt(:,1,2))/nx

! now initialize the velocity field. Note that space has been allocated by CREATE_OPERATORS
u=>Info(1)%ElliInfo%qe(1:nx+1,1:ny,1:nz,1);
v=>Info(1)%ElliInfo%qe(1:nx,1:ny+1,1:nz,2);
CALL RANDOM_NUMBER(u)
CALL RANDOM_NUMBER(v)

! here we link the fields in the second Info to the ones defined in the first. 
DEALLOCATE(Info(2)%ElliInfo%qe) ! first we deallocate the space (to avoid a memory leak)
Info(2)%ElliInfo%qe=>Info(1)%ElliInfo%qe ! we link them. Now anything done by 2 is seen by 1 and viceversa

! 1.1 Off we go
! The general philosophy of ELLIKIT when it comes to solvers is that solvers are initialized, 
! solutions are calculated, and solvers are shut down.

i=POISSON_SOLVERS_THIN(Info(2)%ElliInfo,INIT) ! here the accelerator is initialized. This process uses
                                              ! a direct solution for the horizontal laplacian, which is
                                              ! calculated now
i=POISSON_SOLVERS(Info(1)%ElliInfo,INIT,solver,OP=CSR_PRECOND) ! the iterative solver is initialized. The last
                                                               ! optional parameter is used to set up the 
                                                               ! preconditioner. 


CALL hodge_thin(Info(2)%ElliInfo,1,.false.,ut,vt) ! Here the accelerator is invoked. 
! Notes 
! 1) It works on the second Info
! 2) The second argument is the number of iteration to apply. Here, only one is taken, corresponding to
! and hydrostatic solution. If the leptic ratio is larger than 0.6, a large value of this parameter will
! result in the exact solution being calculated. Of course, here we are interested in the case of 
! small leptic ratios.
! 3) the third argument tells to update u->u-\grad_x P w->w-grad_z P with P being the approximated pressure. If it is set 
! to true, than it calculates the vertical velocity as w=-int (du/dx)dz. In the former case the divergence of the field will
! not be zero, the the latter it will. Of course here, we want the first option.
! 4) ut and vt are the values that we want at the boundaries.





ITERATE=.true.;i=1 ! start the iterations
DO WHILE(ITERATE.AND.i<10000) ! don't do more that 10000 iterations
PRINT*,'iteration ',i
i=i+1

! NOTE: here the projection step with the iterative solver is spelled out rather in detail. In fact
! there is a routine (with and interface similar to hodge_thin) that does the same.

! here we set up the NM conditions for the pressure field.
Info(1)%ElliInfo%bc(1)%p(1:ny,1:nz)=u(1,:,:)-ut(:,:,1) ! u*-utarget
Info(1)%ElliInfo%bc(2)%p(1:ny,1:nz)=u(nx+1,:,:)-ut(:,:,2)
Info(1)%ElliInfo%bc(3)%p(1:nx,1:nz)=v(:,1,:)-vt(:,:,1)
Info(1)%ElliInfo%bc(4)%p(1:nx,1:nz)=v(:,ny+1,:)-vt(:,:,2)
!!$


PRINT*,'calculates divergence'


CALL take_div(Info(1)%ElliInfo,.true.) ! Calculates divergence and apply Boundary conditions.



div=>Info(1)%ElliInfo%div ! short hand
PRINT*,'divergence stat'
PRINT*,'max divergence before = ',maxval(abs(div))

PRINT*,'rms divergence before = ',sqrt(SUM(div**2)/(nx*ny*nz))

PRINT*,'sum of div before =',SUM(div)
!IF(solver==DCG_SOLVER) j=POISSON_SOLVERS(Info(1)%ElliInfo,INIT,solver) ! need to reinit solver at each iteration for this lousy one
Info(1)%ElliInfo%pold=zero
j=POISSON_SOLVERS(Info(1)%ElliInfo,SOLVE,solver,OP=CSR_PRECOND) ! here we seek the approximate solution (note the SOLVE). 
                                                                ! the parameters of the solver are set so that only a few
                                                                ! iterations are actually performed

! subtract gradient
CALL project(Info(1)%ElliInfo,.true.) ! u->u-grad P, applies BC to P
! invoke the thin solver


PRINT*,'calculates divergence'
10 CALL take_div(Info(1)%ElliInfo,.false.) ! false=>do not reapply BC (they are already set)
div=>Info(1)%ElliInfo%div

PRINT*,'max divergence after = ',maxval(abs(Info(1)%ElliInfo%div))

PRINT*,'rms divergence after = ',sqrt(SUM(Info(1)%ElliInfo%div**2)/(nx*ny*nz))



ITERATE=sqrt(SUM(Info(1)%ElliInfo%div**2)/(nx*ny*nz))>1e-10 ! decide weather to iterate again based on residual value of divergence
IF(ITERATE.AND.MOD(i,how_many_steps)==0) CALL hodge_thin(Info(2)%ElliInfo,1,.false.,ut,vt) ! every how_many_steps apply the accelerator
ENDDO ! end of iteration
i=POISSON_SOLVERS_THIN(Info(2)%ElliInfo,FINISH) !wrap up accelerator 
i=POISSON_SOLVERS(Info(1)%ElliInfo,FINISH,solver,OP=CSR_PRECOND) ! wrap up main solver
!!$PRINT*,'parison_mopup'
!!$open(10,file='u.dat',form='unformatted')
!!$write(10) Info(1)%ElliInfo%qe(1:nx+1,1:ny,1:nz,1)
!!$close(10)
!!$open(10,file='w.dat',form='unformatted')
!!$write(10) Info(1)%ElliInfo%qe(1:nx,1:ny+1,1:nz,2)
!!$close(10)
!!$open(10,file='p.dat',form='unformatted')
!!$write(10) Info(1)%ElliInfo%p(1:nx,1:ny,1:nz)
!!$close(10)



NULLIFY(Info(2)%ElliInfo%qe);ALLOCATE(Info(2)%ElliInfo%qe(1,1,1,1)) ! necessary to avoid troubles
DO i=1,2
error=FREEUP_OPERATORS(Info(i)%ElliInfo) ! freeup all the space allocated for the definition of the operators. 



 
DEALLOCATE(Info(i)%ElliInfo)
DEALLOCATE(Info(i)%topo%h,Info(i)%topo%hx)
DEALLOCATE(Info(i)%topo)
END DO
DEALLOCATE(ut,vt)


END ! end of program

